//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : LowPassFilter.cpp
//----------------------------------------------------------------------------

#include "LowPassFilter.h"
#include "assertplus.h"

#ifdef _WINDOWS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#endif

/// constructor
LowPassFilter::LowPassFilter(int Nx1,
                             int Nx2,
                             int half_width,
                             const Array2D<SHORT> &wet)
: m_Nx1(Nx1)
, m_Nx2(Nx2)
, m_halfwidth(half_width)
, m_wet(wet)
, m_wtemp(Nx1 + 3, Nx2 + 3)
{
}

/// destructor
LowPassFilter::~LowPassFilter()
{
}

/// Low-pass filter of effective width 1 + 2 * halfwidth
/** (c) 2006 Albatros Flow Research, Zwolle, The Netherlands, Programmed by: Gert Klopman */
void LowPassFilter::filter(Array2D<REAL> &w_filt)
{
    // this is done inside the WavesComputer::compute() parallelization:
    // filter entire field (parallel for)

    ASSERT(m_halfwidth >= 1, "halfwidth must be larger than 0 but is " << m_halfwidth);
    ASSERT((m_halfwidth * 2) + 1 <= m_Nx1, "halfwidth must be less than (" << m_Nx1 << "-1)/2 but is " << m_halfwidth);
    ASSERT((m_halfwidth * 2) + 1 <= m_Nx2, "halfwidth must be less than (" << m_Nx2 << "-1)/2 but is " << m_halfwidth);

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // filter in x1-direction
    // (m_Nx1 + 2)*(m_Nx2 + 2) assignments
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            m_wtemp[i1][i2] = w_filt[i1][i2];
        }
    }
    // for ---------------------------------------------------------------------

    // m_halfwidth*m_Nx1*m_Nx2 computations
    for (int j = 1; j <= m_halfwidth; ++j)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            REAL w_filt_f = 1.0 / 12.0 * (2 * w_filt[1        ][i2] - 5 * w_filt[2        ][i2] +
                                          4 * w_filt[3        ][i2] -     w_filt[4        ][i2]);
            REAL w_filt_l = 1.0 / 12.0 * (2 * w_filt[m_Nx1    ][i2] - 5 * w_filt[m_Nx1 - 1][i2] +
                                          4 * w_filt[m_Nx1 - 2][i2] -     w_filt[m_Nx1 - 3][i2]);
            REAL w_copy1 = w_filt[1][i2];
            for (int i1 = 2; i1 <= m_Nx1 - 1; ++i1)
            {
                REAL w_copy2 = w_filt[i1][i2];
                w_filt[i1][i2] = -(w_copy1 + w_filt[i1 + 1][i2]) / 4 + w_copy2 / 2;
                w_copy1 = w_copy2;
            }
            w_filt[1    ][i2] = w_filt_f;
            w_filt[m_Nx1][i2] = w_filt_l;
        }
    }
    // for ---------------------------------------------------------------------

    // m_Nx1*m_Nx2 computations
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            m_wtemp[i1][i2] -= m_wet[i1][i2] * w_filt[i1][i2];
        }
    }

    // filter in x2-direction
    // (m_Nx1 + 2)*(m_Nx2 + 2) assignments
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            w_filt[i1][i2] = m_wtemp[i1][i2];
        }
    }
    // for ---------------------------------------------------------------------

    // m_halfwidth*m_Nx1*m_Nx2 computations
    for (int j = 1; j <= m_halfwidth; ++j)
    {
        for (int i1 = 1; i1 <= m_Nx1; ++i1)
        {
            REAL w_filt_f = 1.0 / 12.0 * (2 * w_filt[i1][1        ] - 5 * w_filt[i1][2        ] +
                                          4 * w_filt[i1][3        ] -     w_filt[i1][4        ]);
            REAL w_filt_l = 1.0 / 12.0 * (2 * w_filt[i1][m_Nx2    ] - 5 * w_filt[i1][m_Nx2 - 1] +
                                          4 * w_filt[i1][m_Nx2 - 2] -     w_filt[i1][m_Nx2 - 3]);
            REAL w_copy1 = w_filt[i1][1];
            for (int i2 = 2; i2 <= m_Nx2 - 1; ++i2)
            {
                REAL w_copy2 = w_filt[i1][i2];
                w_filt[i1][i2] = -(w_copy1 + w_filt[i1][i2 + 1]) / 4 + w_copy2 / 2;
                w_copy1 = w_copy2;
            }
            w_filt[i1][1 ] = w_filt_f;
            w_filt[i1][m_Nx2] = w_filt_l;
        }
    }
    // for ---------------------------------------------------------------------

    // m_Nx1*m_Nx2 computations
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            w_filt[i1][i2] = m_wtemp[i1][i2] - m_wet[i1][i2] * w_filt[i1][i2];
        }
    }
    // for ---------------------------------------------------------------------
}

